
# bitcount taken from mpmath by F. Johansson
from bisect import bisect
powers = [1<<_ for _ in range(300)]

def bitcount(n):
    """Calculate bit size of the nonnegative integer n."""
    bc = bisect(powers, n)
    if bc != 300:
        return bc
    bc = int(math.log(n, 2)) - 4
    return bc + bctable[n>>bc]

# Used to avoid slow function calls as far as possible
bctable = map(bitcount, range(1024))


# part of sympycore/arithmetic/number_theory.py
# Copyright (c) 2007-2008 Pearu Peterson
#
# multinomial_coefficients modified by Mario Pernici

def binomial_coefficients(n):
    """Return a dictionary containing pairs {(k1,k2) : C_kn} where
    C_kn are binomial coefficients and n=k1+k2."""
    d = {(0, n):1, (n, 0):1}
    a = 1
    for k in xrange(1, n//2+1):
        a = (a * (n-k+1))//k
        d[k, n-k] = d[n-k, k] = a
    return d

def binomial_coefficients_list(n):
    d = [1] * (n+1)
    a = 1
    for k in xrange(1, n//2+1):
        a = (a * (n-k+1))//k
        d[k] = d[n-k] = a
    return d

def code_t(t,bits_exp):
  m = len(t)
  expv = 0
  for i,n in enumerate(t):
    expv += n<<(i*bits_exp)
  return expv


def multinomial_coefficients(m, n, _tuple=tuple, _zip=zip):
    """Return a dictionary containing pairs ``{(k1,k2,..,km) : C_kn}``
    where ``C_kn`` are multinomial coefficients such that
    ``n=k1+k2+..+km``.

    For example:

    >>> print multinomial_coefficients(2,5)
    {(3, 2): 10, (1, 4): 5, (2, 3): 10, (5, 0): 1, (0, 5): 1, (4, 1): 5}

    The algorithm is based on the following result:
    
       Consider a polynomial and it's ``m``-th exponent::
       
         P(x) = sum_{i=0}^m p_i x^k
         P(x)^n = sum_{k=0}^{m n} a(n,k) x^k

       The coefficients ``a(n,k)`` can be computed using the
       J.C.P. Miller Pure Recurrence [see D.E.Knuth, Seminumerical
       Algorithms, The art of Computer Programming v.2, Addison
       Wesley, Reading, 1981;]::
       
         a(n,k) = 1/(k p_0) sum_{i=1}^m p_i ((n+1)i-k) a(n,k-i),

       where ``a(n,0) = p_0^n``.
    """
    m = int(m)
    n = int(n)
    if m==2:
        return binomial_coefficients(n)
    bits_exp = bitcount(n)
    mask_exp = (1<<bits_exp)-1
    masks = [1 << (bits_exp*i) for i in range(m)]
    symbols = [(0,)*i + (1,) + (0,)*(m-i-1) for i in range(m)]
    s0 = symbols[0]
    p0 = [_tuple(aa-bb for aa,bb in _zip(s,s0)) for s in symbols]
    p0 = [code_t(t,bits_exp) for t in p0]
    r = {_tuple(aa*n for aa in s0):1}
    r_get = r.get
    r_update = r.update
    l = [0] * (n*(m-1)+1)
    l[0] = [(code_t(_tuple(aa*n for aa in s0),bits_exp),1)]
    #print 'DB0', l[0]
    for k in xrange(1, n*(m-1)+1):
        #print 'DB5 l[%d]=%s' %(k-1,l[k-1])
        d = {}
        d_get = d.get
        for i in xrange(1, min(m,k+1)):
            nn = (n+1)*i-k
            if not nn:
                continue
            t = p0[i]
            #print 'DB1 t=',t
            for t2, c2 in l[k-i]:
                #print 'DB2 t2=',t2
                #tt = _tuple([aa+bb for aa,bb in _zip(t2,t)])
                tt = t2 + t
                #print 'DB3 tt=',tt

                cc = nn * c2
                #print 'DB6', bitcount(tt), bitcount(cc)
                d[tt] = d_get(tt,0) + cc
        b = []
        for t,c in d.iteritems():
          if not c:
            continue
          c //= k
          if (t&mask_exp) > 0:
            b.append((t,c))
          a = [0]*m
          i = 0
          while t:
            exp = t&mask_exp
            t >>= bits_exp
            a[i] = exp
            i += 1
          #tdec = tuple(a)
          r[_tuple(a)] = c
        l[k] = b
    return r

def newton_method_sizes(N):
    """taken from Sage misc.py
    """
    N = int(N)
    if N < 1:
        raise ValueError, "N (=%s) must be a positive integer" % N

    output = []
    while N > 1:
        output.append(N)
        N = (N + 1) >> 1

    #output.append(1)
    output.reverse()
    return output

